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I. INTRODUCTION 

It is well known that the asymptotic behavior of the two-point function for large Euclidean 
times is determined by the lowest eigenvalue of the Hamiltonian in a given channel. In case 
of stable particles, this property allows one to determine their masses. The case of the 
excited states is different. Here, the two-point function yields the spectrum of the so-called 
"scattering states." The relation to the energy and width of the resonance states is not 
direct, since a resonance, in general, can not be associated with an isolated energy level of 
a Hamiltonian. Up to now, several alternative methods have been used to determine these 
quantities from the lattice Monte Carlo (MC) simulations. These are: 



At present, Liischer's approach [iH^ is widely used to deal with the resonances in 
lattice QCD, obtained in simulations with sufficiently low quark masses. In brief, 
the procedure consists in determining first the phase shift by studying the volume 
dependence of the energy spectrum on the lattice. Then, continuing the S'-matrix 
into the complex plane (e.g., by using the effective-range expansion whenever pos- 
sible), one attempts to determine the position of thepoles on the second Riemann 
sheet. This procedure is described in detail in Ref. (5|, where the generalization to 
the case of the resonance matrix elements (in 1-1-1 dimensions) is also considered. A 
shortcut is provided by using Breit-Wigner type parameterization for the scattering 
phase and determining its parameters (energy and width) from the lattice data (see, 
e.g. Ref. pI and Refs. [7|, |8[, where the method has been applied in the case of the 
a- and p- mesons, respectively). The Liischer's approach has been also generalized for 
the moving frames |9(]. 

Recently, the spectral functions in QCD have been reconstructed by using the maximal 
entropy method (see, e.g. 10l-ll2|). This method, as well as Liischer's approach, has 



in principle the capacity to address the problem of the extraction of the resonance 
energy and width from the Euclidean MC simulations on the lattice. 

iii) The Euclidean correlators have been parameterized in terms of the energy and width 
of an isolated resonance state, in order to subsequently determine these quantities 
from the fit to the lattice data 13|. In that paper, the method has been applied to 



study the glueball decay, 
iv) In certain cases, the decay width of an excited state can be evaluated by calculating 



decay amplitudes on the lattice (see, e.g. |14|. Il5l|). 



Recently, there has been a substantial activity in the determination of the excited 
meson and baryon spectrum by using generalized eigenvalue equations 16l-l23| . Despite 



spectacular progress achieved in the field, it should be stressed once again that a 
resonance state can not be uniquely associated with a particular energy level. To 
a certain extent, excited states and scattering states can be distinguished, e.g., by 
studying the volume dependence of the spectral density |2J, |25| . This method, however 
works for narrow resonances only |26| . 



In this paper, we combine some of the above ideas and propose a systematic method 
to extract resonance pole positions from lattice data. In its present form, our approach is 
applicable to the systems with a low-lying, well-isolated, narrow resonance in the spectrum 
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FIG. 1: Schematic representation of the analytic structure of the momentum-space two-point 
function in the rest frame p^ = (f^, 0) as a function of the complex variable s = (iuj)'^, for two 
different values of the box size L. The crosses on the real axis denote the poles, and the shaded 
blobs, which are located symmetrically to the real axis, mark the location of the resonance poles 
on the second Riemann sheet emerging in the infinite volume limit. At any finite value of L, 
however, the two-point function is meromorphic and the second sheet does not arise. If the two- 
point function is measured for Euclidean times < t, the "energy resolution" of such a measurement 
is approximately equal to t~^, as indicated on the right panel. 

(for example, the p or the A (1232)). First, we have tested our method using synthetic 
input data, represented by the Euclidean propagator of the A, calculated in the low energy 
effective field theory at one loop. A further test has been carried out in a 1+1 dimensional 
model of two coupled Ising spins, where the resonance parameters have been determined in 
the past utilizing Liischer's approach |27|, |28|. In both cases, we find that the method is 



capable to extract the pole position of the resonance. 

The outline of the paper is as follows. In section [TT] we discuss the foundations of the 
method. The general representation of the two-point function in the presence of a low-lying 
isolated resonance is discussed in section IIIII In section HV] we consider the procedure of 
the data fitting and the determination of the pole position by using synthetic data. A short 
review of the 1+1 dimensional Ising model is given in section |Vl The extraction of the 
resonance pole in this model is considered in section I VII Finally, section IVIII contains our 
conclusions. Some technicalities are relegated to the appendices. 



II. KALLEN-LEHMANN REPRESENTATION 



In the beginning of this section, we present a qualitative reasoning to justify our method. 
We start by mentioning that lattice QCD simulations are always carried out on lattices with 
a finite (Euclidean) time and spatial extension. Below, we consider lattices of the size TxL'^, 
where T and L denote the box size in time and space, respectively, and d is the number of 
spatial dimensions. For not so large L (however, large enough to suppress the polarization 
effects in stable particles), the energy levels are well separated and can in principle be 
extracted from the asymptotic behavior of the two-point function at large Euclidean times 
t ^ AE, where AE denotes the average level spacing. The Fourier-transform of the two- 
point function is a meromorphic function in the complex s-plane (see Fig. [T^). Consequently, 
the second Riemann sheet, as well as the poles on it (corresponding to the resonances), do not 



appear at any finite L. The information about tliese poles stays encoded in tlie dependence 
of the spectrum on the spatial box size L and can be extracted (in several consecutive steps, 
as explained in the introduction) by using Liischer's approach. 

In case an almost stable state is present, such a complicated procedure might seem 
superfluous. Intuitively, it is clear that, if the decay width F is small, the resonance will 
behave pretty much the same way as a stable particle and will determine the t-dependence 
of the two-point function in a large interval (however, not for the asymptotically large times 
t ^ r~^, when the resonance already has decayed). Consider, for example, a lattice with 
the large spatial size L, and calculate the two-point function at finite t (see Fig. [Jd). If the 
"resolution" ~ t~^ is larger than the distance between energy levels, the two-point function 
is given by a sum of (many) exponentials with the spectral weight suppressed by a factor 
L~'^, so that, effectively, the spectral sum transforms into an integral over the energies. This 
is equivalent to the emergence of the cut that connects the physical sheet to the second sheet. 
We expect that, in this case, one may find an alternative representation of the two-point 
function in terms of quasiparticle degrees of freedom corresponding to the pole(s) on the 
second Riemann sheet (i.e. the resonance decay and width) plus a small background which 
can be described by a few parameters. Namely, we expect that there exists a window in 
the variable t, where such a description will be more effective than the multi-exponential 
representation through the energies and spectral weights. 

Hence, the original problem is reduced to finding a universal, model-independent param- 
eterization of the two-point function in the presence of a narrow resonance, which will allow 
one to determine the energy and the width of the latter by performing a fit to the lattice 
data. This is analogous to the exponential parameterization, which allows one to determine 
the mass of a stable particle by fitting at asymptotically large times. Note that we shall 
try to avoid approximations in the spectral function (e.g. the narrow width approximation 



used in Ref. 13|), which are a potential source of a systematic error. The contribution of 
the background, albeit small, will be taken into account in a systematic manner. 

The model-independent parameterization, which was mentioned above, can be obtained 
directly by using the Kallen-Lehmann representation for the two-point function. Below, we 
give a brief derivation of this representation in a finite box. Further, we consider the limits 
T — )■ oo and L — > oo in detail, in order to quantify the qualitative arguments given in the 
beginning of this section. 

In the derivation of the Kallen-Lehmann representation in a finite volume, we mainly 



follow the steps given in Ref. [10|. Let 4>{x) = 0(t, x) be any local field, carrying the 
resonance quantum numbers. The two-point function of this field in the Euclidean space is 
defined as 

D(t, x) = {N, [0(t, x)0t(o)] ) = ^ Tr {Nt [0(t, x)0t(o)] e"^^) , 

iV,[0(t,x)0t(o)] = ^(t)0(t,x)0t(o)±^(-t)0t(o)0(t,x), 

ZiT) = Tr (e-^^) , (1) 

where the upper (lower) sign corresponds to bosons (fermions). The Fourier transform of 
this expression takes the form 

D{t, x) = ^ 5^ 5^ D{zu, k) e— ^-^'^ , (2) 

U! k 



where to = ^ut {co = ^ (2raT + 1)), with ut G Z, are Matsubara frequencies in case of 
bosons (fermions), and k = ^ n with n e Z"^. 

Using a complete set of Hamihonian eigenvectors H\a) = Ea\a) to calculate the trace in 
Eq. ([1]), one gets 

^(^^>k) = -^ EE<p.-p. '/- I'+TJ {a\mm{f^mo)\a) , (3) 

where (5^ denotes the periodic Kronecker 6 in d dimensions. This relation can be rewritten 



as 



-r^^—Aiu',k), (4) 

where the spectral function is given by 

A{u;',k) = -^ EE(^"''"^T^"''''^)'^(^'-^/^ + ^")<P.-P. («l0(O)|/3)(/3|0t(O)|a).(5) 

By applying discrete symmetries, it can be shown that the spectral function obeys the 



following properties [10 1 



A(w',k)>0 for a;'>0, A(-w', -k) = ^^(a;', k) = ^^(a;', -k) . (6) 

The dispersion integral can be rewritten as^ 

D(za;,k)= / ^J^A(a;',k). (7) 

Finally, in the limit T — t- oo, only the vacuum state a = 0, i^a = contributes, and the 
spectral function in Eq. ([7]) is given by 

^hm A(a;',k) = L' E^(^' "^/^Xp. mm\P)? ■ (8) 

Now, let the variable cu^ be outside the narrow strip along the negative real axis. Then, the 
function 2uj'{uj' + o;^)"^ for all < cj' < oo is uniformly bound from above by some large 
constant B. For fixed B, performing the limit L — )■ oo and applying the regular summation 
theorem [29|, it is seen^ that the quantity D{iuj, k) approaches D°°(iuj, k), with the pertinent 
spectral function given by 

A-^ico', k) =^ Sico' - E„) {2iiY5\]^ - p^) I (O|0(O) 1/3) p , (9) 



^ For illustrative purpose, below we display the bosonic case only. The fcrmionic case can be treated 

similarly. 
^ The regular summation theorem implies that the matrix elements (O|0(O)|/3) are continuous functions of 

Ep. We examine these matrix elements explicitly in a simple quantum mechanical model in Appendix VK\ 

and show that in this case the above requirement is indeed fulfilled. 



where ^ « stands for the sum (integral) over the continuous spectrum wave functions. More 
precisely, for uj^ outside the strip 

/•oo 

\D{tuj,'k)-D'^{icu,'k)\<B du'\A{u','k)-A°°{u',k)\, (10) 

where Umm is determined by the invariant mass of the lowest-mass state. Further, according 
to the regular summation theorem, the difference A{uq, A, k) — A°°{uq, A, k), where 



/.wo+A/2 /.wo+A/2 

A{uJo,A,k)= duj'A{io',k), A°°(wo,A,k)= / dio'A^{uj',k), (11) 

Jujo-A/2 JuJn-A/2 



for any uq > Wmin and A > converges faster than any power of L as L — > oo (for a 
discussion, see Appendix |A]) . In other words, in this limit the two-point function converges 
to its infinite-volume counterpart everywhere in the complex plane except the narrow strip 
along the cut (for a related discussion, see also Ref. |30|). This statement is a mathematical 
formulation for the intuitive picture of "poles merging into the cut," which is shown in 
Fig. [Do. To further illustrate this, an example of a function which is meromorphic at a finite 
L and develops a cut and a pole on the second Riemann sheet in the limit L — ;■ oo, is given 
in Appendix [Bl 

Moreover, from Eq. ([2]) one finds 

D{t, k) = -J2 ^-""-'m^, k) = / du' _ +,^ A{u', k) 

POO 

-^ / duj'e-^''A{uj',k). (12) 

The last line is obtained in the limit T — t- oo. Together with the expression ([8]) for the 
spectral density, we recover the representation for D{t,k.) as a sum over exponentials. In 
the case t~^ is much larger than the distance between different energy levels (this can be 
achieved, e.g., by holding t fixed and increasing L), many exponentials contribute to D{t, k) 
and the sum over the energy eigenvalues can be replaced through the integral. In this case, 
A{uj', k) is replaced by A°°{uj', k). 

To summarize, the behavior of the two-point function can be studied in different regimes. 
For asymptotically large t and moderately large L, only the few lowest, well-separated energy 
levels contribute. This situation is well described by a sum of a few exponential terms. In 
difference to this, in the regime with asymptotically large L and moderately large t there 
are many terms with nearly the same energies that contribute to the multi-exponential 
representation. The sum over the discrete energy spectrum effectively transforms into an 
integral. If, in addition, a low-lying well separated resonance emerges, we expect that 
the spectral integral can be efficiently parameterized in terms of the resonance parameters 
instead of the stable energy levels. 

III. TWO-POINT FUNCTION AT FINITE t 

As discussed in the previous section, it is possible to perform the infinite- volume limit L — )■ 
oo in the two-point function, keeping the Euclidean time t fixed. The spectral representation 
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is given by Eq. ([7]) with the spectral density given by Eq. (Q. Note that the spectral density 
vanishes for < oj' < Wm\n, and hence the integration in Eq. ([7]) in fact is performed from 
^'^ = ^min to infinity. 

For simplicity, we work in the center-of-mass (CM) frame k = and denote D{iu}, 0) = 
D{iuj), A°°{uj',0) = A{uj'). The spectral representation then takes the form 

DM=f i^^M. (13) 



min 



In the vicinity of the elastic threshold, A(uj') ~ (tu' — LOmmY'^^^'^, where / stands for the orbital 
angular momentum^. 

Assume now that an isolated low-lying resonance emerges. This is equivalent to the 
statement that the function A{u') takes the form 

where the singularities of the function Q{uj') lie far enough from the threshold so that the 
Taylor expansion of this function converges in the part of the complex plane that includes 
the resonance poles at u' = u^ and u' = u*^. The energy and the width of the resonance 
is determined by ur in the standard manner. Note that these poles come in complex- 
conjugated pairs'^. 

Using Eq. (TT^ . it can be easily shown that 

/•oo /"OO JTpTpl+l/2„-Et 

D{t) = l duj'e-''A{J)=e--' j^ {E - En){E - E^) ^^^ ^ """"'-^ ' ^^^^ 

where E^i = ur — Umm = Eq — iT/2. As mentioned before, it is assumed that the Taylor 
expansion Q{E+Umm) = YlkLo QkE'^ converges in the part of a complex region which includes 
the resonance poles. 

From the above expression we get 

°° /-oo j-rp-rprn+l/2 -Et 

D{t) = e--"^^q,F(>^+^\t,En), F(-\t,En) = j^ ^^ _ ^^^^ ^^^^^ . (16) 

In particular, for m = 0, 1, we find 

F^'\t,ER) = -l,lmx{t,En) 

F^'\t,En) = Rex(t,i?if) - ^ Imx(t, i^^), (17) 



where 



x{t,E^) = J^ ^_^^ (18) 



^ This statement is valid in 3+1 dimensions. In 1+1 dimensions, one has to substitute / = in all formulae. 
^ A{uj') is the discontinuity of a function which is analytic in the cut complex plane and obeys Schwarz 

reflection principle. Hence the poles in this function (which emerge on the second Riemann sheet), always 

come in pairs. This is the justification for the ansatz (|14p. 




FIG. 2: Schematic representation of the effective mass for a stable particle (dashed line) and for a 
resonance (solid line). 



and the following representation in form of an infinite series is useful in a wide range of the 
variable t 



x{t,Ef 



-n 



-E 



Re 



-ErI 



oo 

E 

fc=0 



-2ERtf 



k+l 



{2k 



(19) 



Further, the functions F*-™-* with m > 2 can be recursively expressed through F*-°'^\ The 
general representation of the two-point function follows straightforwardly from the above 
relations 



y- OO 

D{t) = e---"* I coF(o)(t,ER) + c,F^'\t,Ej,) + ^ 

'^ fc=0 



Xk \ 

fl+k+3/2 J ' 



(20) 



where co,i and Xk are expressed through the Taylor coefficients qk as well as EQ,r. 

Eq. fl20|) represents our central result. It provides a universal parameterization of the 
Euclidean two-point function in the presence of a low-lying isolated resonance described by 
two parameters Eq, F. The couplings x^ are associated with the non-resonant background. In 
particular, it encodes the contribution of the threshold which lies below the resonance energy. 
This means that if t is taken too large, the background dominates and the information about 
-Eq) r is erased. We assume, however, that in the presence of a narrow resonance, there exists 
a sufficiently wide window in t, where the background is small and Eq, F can be determined 
from the fit of the measured D{t) to the representation fl2U]) . We require that, in this window, 
adding the background parameterized by the constants Xk should lead to small corrections 
in Eo,r and the fit should remain stable against the increase of the number of independent 

Xk. 

The physical meaning of our method can be easily illustrated by Fig. [2l In this figure, 
the effective mass of a system in the presence of a stable state/resonance is schematically 
depicted. If there is a stable particle, the plateau in the effective mass sets in almost 
immediately. However, if a resonance instead of a stable particle is present, there exists a 
wide window in t, where many excited states contribute and the effective mass decreases 
slowly until it reaches the asymptotic value. Our method roughly corresponds to fitting 
D{t) within this interval by the representation given in Eq. fl20|) . and the decay width is 
determined by the rate of the decrease of the effective mass. Note that a similar picture 
was obtained in Ref . 31[ . In that paper, the theory of two coupled scalar fields, where the 
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heavier field can decay into a couple of light scalars, was considered. In particular, it has 
been shown that the effective mass of the heavy scalar, calculated on the lattice, exhibits 
the same behavior when its mass is below/above the two-particle threshold (see Fig. 1 of 
that article) . The two point functions of the excited mesons in QCD also exhibit a similar 
behavior 16| . 



It is instructive to compare the parameterization fl2UI) to the pertinent formula obtained 



in Ref. 13|. A technical difference consists in the absence of the threshold factor [oj' — 
'^min)^^^ that results in a different parameterization of the background. The important 



difference is, however, that in Ref. [13| a Breit-Wigner type parameterization for the spectral 
function is originally assumed. Since this ansatz did not fit the lattice data well, an ad hoc 
energy-dependence in the decay width has been introduced, and the functional form of 
this dependence has been determined by using a trial-and-error method. In our approach, 
the representation of -D(t), given in Eq. fl2U]) is completely general and is based on the 
sole assumption that a well-isolated resonance emerges at low energies. The background is 
parameterized by the constants Xk in a systematic manner. 

IV. THE FIT 

We first test our method by using synthetic data. Consider, for instance, the propag ator 



of the A-resonance evaluated in the Small Scale Expansion (SSE)^ at one loop [33|, |3J]. In 
the Minkowski space, this propagator is given by 

5a (p) = -5f (p)n3/2 + spin-l/2, 

^^ ^^^ ^ mA(l + S2(p2))-^(l-Si(j92))' ^^^^ 

where rh/^ denotes the mass of the A-resonance in the chiral limit, 11^/^ stands for the 
projector onto the spin-3/2 state and the spin-1/2 part does not have a pole in the low- 
energy region (for a general proof of this statement, see |35|). Further, the invariant functions 
^1,2 (p^) at order p^ in the chiral expansion are given by 

Si(^) = % {W,{s) - W,{s)) , T.,{s) = ^!I^ W^is) , (22) 

where ca and F denote the vrA^A coupling constant and the pion decay constant in the chiral 
limit, respectively, mN,m^,MT^ are the nucleon, A and pion masses, respectively, and the 



^ The SSE is a phenomenological extension of Chiral Perturbation Theory in which the A-nucleon mass 
spUtting is counted as an additional small parameter. This quantity, however, does not vanish in the 



chiral limit. The framework of the SSE is laid in detail in Ref. 



32| 



invariant functions 1^2,3 are given by J33|, |34 



_ -s^,n%-M; ^ _Mi / _^ M3 1 






sTA s-ml + M^ f M^ \ ^ s + Ml-m\ + ^/\ 



°^ ^ 167rs 327r2s V ^^ / 327r2s ^ + M2 - m^, - ^/X ' 

A = (s - (m^ + M^)%s - {rriN - M^f) . (23) 

The trace of Sj^ (p) obeys the dispersion relation 

Tr S'/'h) - 4m^(l + S2(p^)) _ r ds' 

where the expression for the discontinuity can be directly read off from Eqs. fl2Tl) - fl23l) . 

In the calculations we have used the following values of the parameters: ttln = 940 MeV, 
M^ = 140 MeV, rriA = niA = 1232 MeV, F = F^ = 92.4 MeV and ca = 1.5 (this value 
leads to the width F = 124 MeV in a O(e^) calculation at p^ = m\). It is easy to check that 
the propagator has a pole at rriR = 1212 MeV and F = 76 MeV (note the large shift in the 
quantity F as compared to its value obtained at p'^ = m\ that presumably is an artefact of 
a 0{e^) approximation). 

Next, we wish to investigate whether it is possible to recover this result by applying our 
method. To this end, we analytically continue Eq. (!24|) into Euclidean space and perform the 
Fourier transform with respect to the fourth component of the momentum. The resulting 
values are treated as synthetic data. We choose the interval 1.7M~^ < t < AM~^ and 
perform a least squares fit of these data to the formula (120]) (the data points are assumed 
to be distributed equidistantly in this interval). 

In the fit, we cut the sum in Eq. fl20j) at some value k^ax- The fit of the 7 data points 
with kmax = yields mn = 1213 MeV and F = 74 MeV that is already close to the exact 
values. The procedure converges rapidly. At the accuracy of the digits displayed, the exact 
result is obtained for /c^ax = 2. Adding more terms, it is possible to improve the agreement 
with the exact result up to very many decimal digits. 

To summarize, using synthetic data, we have demonstrated that our method is capable 
to reconstruct the exact position of a pole in a complex plane from a limited data sample. 
To perform a similar analysis for real Monte Carlo data is much more challenging. One of 
the main problems that we have encountered there, is related to the instability of the fit 
when kmax increases (this problem already arises for relatively small kmax = 3 or 4). Namely, 
the constants Xk, which describe the background, become very large in magnitude having 
alternating signs and this destabilizes the values of i?o, F extracted from the fit. 

In order to circumvent this problem, we have performed a Bayesian fit to the lattice MC 
data. A detailed description of the Bayesian fit techniques, which is well suited for our 
purposes, can be found, e.g. in Ref. [36|. We shall present a brief summary of the method 
below. The function to be minimized in the standard least squares fit is given by 

X" = y2{D{U, Eo, F, ci, C2, Xk) - D{U)f , (25) 
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where D{ti) are data corresponding to the points tj. In Eq. (125|) it is imphcitly assumed 
that the MC errors in the data D(ti) do not vary much with tj. Note that the above form 
still does not include our prior knowledge about Xk- The assumption about the smoothness 
of the function Q{u') in Eq. ( 1T^ implies that Xk should be of "natural size" excluding the 
scenario where the Xk become large with alternating signs. 

In order to implement this prior knowledge into the fitting procedure, in analogy with 



Ref. [36|, we define the augmented x 

y2 = y2 . y2 2 ^ J_ ST^ 2 (2Q) 

Aaug A I Aprior ' Aprior 02 / ^ k ' \ / 

k=0 

where 5* is some scale that ensures that all Xk stay in the "natural" range. 

We determine the quantity S by using the trial-and-error method. If S is too large, the 
introduction of xLg does not cure the problem with the convergence. This sets the upper 
limit on the value of S. The lower limit for S is set by the requirement that the results 
obtained with standard x^ and xLg agree for low k^ax = 1,2. In addition, within this range, 
the final result of the fit for Eq, T should not depend on S. 

In the 1+1 dimensional model with two Ising spins discussed in the next section, we 
have performed fits using xLg- Below we show that this technique allows one to extract the 
precise values of i^o, T from the lattice MC data in this model. 

V. 1 + 1 DIMENSIONAL MODEL WITH TWO COUPLED ISING SPINS 

In this section we apply our method to the extraction of the resonance pole position 
to a 1+1 dimensional model of two coupled Ising spins. This model has been treated 



in Refs. |27l . |28| using Liischer's approach. In particular, it has been shown that a narrow 
resonance emerges in the system, whose parameters can be extracted in a systematic manner. 
The action of the model is given by 

S = -K^ ^ (pz(pz+f,- Hr, J^ VzVz+f, + - ^ Vz(pz{(pz-f, + (t)z+il), (27) 

zeA,fi=l,2 z<=A,fi=l,2 zeA,fi=l,2 

where (pzjVz = =t 1 are two Ising spins which interact with each other through the Yukawa- 
type coupling gri(f)(f). The sum 2; G A, where z = (x, t), runs over all lattice points and fi 
denotes the unit vector along the spatial axis. The couplings k,^, k^ > are chosen so that 
the masses of and rj are m^ ~ 0.19 and m^ ~ 0.5 (in lattice units). Note that, ii g ^ 0, 
the rj decays into 2(f), so m^ corresponds to the resonance energy in this case. 



The model has been analyzed in detail in Refs. [27|, |28|. We give only a short summary 
of this analysis here. In particular, it has been argued that in the theory described by the 
Lagrangian fl27|) no second-order phase transition occurs and thus the continuum limit can 
not be performed. In other words, all results obtained here refer to the effective theory with 
an ultraviolet cutoff. 

The energy spectrum is determined by solving the generalized eigenvalue problem. The 



operator basis is defined as in Refs. |27l . |28 



Oiit) = jY,V.,t, C?.W = ^$^0x,t0,,te-^^(^-^), P. = ^^^V^' (28) 
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FIG. 3: Checking the convergence of m^ (left panel) and T^j (right panel) against the variation of 
^max and ^max) for L = 60. A similar behavior is observed for smaller values of L. 



with J = 2, 3, ■ ■ ■ . The correlator matrix is given by 

a,(t) = {O,{t)O,{0)) - {Om){OM) ■ (29) 

The spectral decomposition of Cij{t) is approximated by the truncated series 



a,{t) = Y,vi 



i%f\-w,t 



(30) 



1=1 



The energy eigenvalues Wi for / = I,--- ,r are determined by diagonalizing the matrix 
Mjt.tn) = C~^'^{to)C{t)C~^''^{to), where to is some fixed time (in the following, as in 
Refs. |27|, |28|, we always use to = !)• The eigenvalue equation takes the form 



M(t,to)u(') = X^^\t,to)u 



(0 



X^'\t,t, 



-Wiit-to) 



l = l---r 



(31) 



where u^'^ form an orthonormal basis. The eigenvectors f*^'^ are given by i;'-'-' 



exp{Wit/2)C^/^{to)u^^l 

The MC simulation is done by using a cluster algorithm 37|. We closely followed the 
procedure described in |27|, |28| and, using the parameter set k^ = 0.3700, k^ = 0.3700, g = 
0.04, have reproduced the L-dependent spectrum calculated in this paper. The resonance 
parameters found in Refs. |27|, l28| are: m^ = 0.5112(3) and F^ = 0.0100(3). It remains to 
be seen, whether the same result can be obtained by using our approach. 



VI. RESULTS 

In order to use our method, one has to calculate the two-point function within a suffi- 
ciently large interval in the Euclidean time t and then fit the result with Eq. (1201) . To this 
end, the correlator Cu(t) has been chosen, see Eq. (^U^. However, as already mentioned 



in [27|, l28|, the simulations become unstable already at t ^ tunst = 5 — 8, depending on 
the value of L chosen. The statistical error in the effective mass of 77 at t > tunst blows 
up, rendering an accurate fit impossible. The use of improved estimators or a substantial 
increase of the number of configurations results only in a moderate improvement of the error 
in Cii(t). 



12 



As described above, the energy spectrum of the system can be determined with high 
accuracy from the correlator matrix Cij{t) at t < tunst by applying the generalized eigenvalue 
method. In addition to the ground state, the approach allows a reliable extraction of higher 
excited levels (up to 4 or 5 levels, depending on L). The physical reason for this is that the 
matrix Cij{t) contains much more information about the system than the single function 
Cii(t). In particular, it contains information about the matrix elements describing the 
transitions between various energy levels. 

So, it is not surprising that using this input in our method helps to reduce the errors 
dramatically and to stabilize the fit. In brief, the procedure can be described as follows: 

1. The energy spectrum Wi and the wave functions f^'^ are accurately determined by 
measuring the matrix Cjj(t) at t < tunst- We choose t^nst = 5 for all L and average all 
ly, fort = to---(tunst-l). 

2. The function Cii(t) is approximated by the multi-exponential function Cii(t) = 
^[^^ 2;/exp(— VF^t), where the zi = |fi P are averaged for all t = to ■■■ (tunst — !)• 
This approximation is used for t > tunst as well. Note that zi, 1^1 encode the in- 
formation about the overlap of t] and 20 states that determines the decay width of a 
resonance. 

3. The expression (120 p is fitted to the Cuit) which is approximated by the multi- 
exponential function. 

The MC simulations were carried out for various lattice sizes in the interval L = 24 — 60, 
while the value T = 100 remained fixed throughout the simulations. We have used bases 
containing 4-6 operators and performed test runs for some (large) values of L by using the 
basis of 8 and 10 operators. In the fit, all data between t = 1 and t = t^ax were used. The 
errors in our results are purely statistical and were estimated by performing 5 independent 
simulations with 10^ configurations each. In addition, we find that the increase of the number 
of operators to 8 or 10 operators does not affect the result within the errors. 

First, the stability of our results was checked, when tmax and fcmax increase. The result 
of this check is displayed in Fig. [3l where the dependence of the real and imaginary parts 
of the resonance pole position on tmax is plotted for different values of fcmax- It is seen that 
for tmax > 10 — 12 both the energy and the width remain almost constant and converge 
rapidly in fcmax already at fcmax = 3, if the Bayesian fit is performed. The similar behavior 
is observed at all values of L. The final result for the resonance pole parameters is always 
given at A;max = 10. 

In order to ensure that, performing the Bayesian fit, a bias is not introduced in the 
extracted values of the resonance parameters, one has to check that there exists a range of 
the scale parameter S where the energy and the width depend weakly on S. The results 
for both quantities at different values of L look qualitatively similar. In Fig. |4] we present 
the plot for the width at L = 60. As seen from Fig. HJ a wide plateau emerges around 
S ~ 10^, where the scale dependence practically disappears while the convergence in /cmax 
still persists. This is the window, where the extraction of the width is finally carried out. 
Increasing S even further, the convergence in /cmax breaks down, and the result can not be 
trusted any longer. 

Finally, since our MC data have been calculated at a finite L, whereas the formula f l20|) 
refers to the limit L — )■ oo, there is an expected residual volume dependence in the parameters 
Eq, r. The Fig. [5] displays this dependence. In particular, it is seen that there is a rather 
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FIG. 4: Dependence of the width on the scale S used in the Bayesian fit at different fcmax- For 
smah values of 5, the result is scale-dependent. For large 5, the result does not converge with 
^max- There exists a plateau around S ~ 10^ where the procedure converges and yields a scale- 
independent result. 

strong variation of the width at small values of L that flattens around L = 48. In the present 
paper we do not attempt to quantitatively describe the finite volume artefacts. This issue 
forms the subject of a separate investigation and we plan to address it in the future. 

From Fig. [5] it is also seen that the effect of the background on the real part of the pole 
position is small, whereas the imaginary part is far more sensitive to it. Namely, the two 
values of Eq, calculated at L = 60 for kmax = and k^ax = 10 differ by ^ 0.2%, whereas the 
same calculation for the width yields F^ = (0.91 ± 0.04) ■ 10-^ and F^ = (1.17 ± 0.05) ■ lO'^, 
respectively. In general, one may conclude that the effect of the background can not be 
neglected. 

The final result for the real and imaginary parts of the pole position (for L = 60) are 



rur, = 2m^ + Eo = 0.5074 ± 0.0004 , 
F^ = F = (1.17 ±0.05)- 10^2 



(32) 



(errors are only statistical). This result can be checked by using the effective-range expansion 
for the scattering phase (cf. with Ref. |27| ) 



— — tanS{p) = a — bp^ 



W = 2jm 



pz 



(33) 



where the parameters a and b are related to m„ and F„ through 



m^ 



2i ml + 



r. 



bm'^ 



(34) 



As one sees from Fig. El our phase shift results are generally in agreement with the results 
of the Ref. |27| . However, since the data are not exactly linear, the question arises, which 
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FIG. 5: The energy and the width of the resonance, extracted from the data at fcmax = 10 (final 
result). For comparison, we display the result at fcmax = (no background) and the result, taken 
from Refs. J27l. |28| (the error quoted in these references corresponds to the thickness of the lines). 



The errors in our calculations are purely statistical. 



interval in the variable p^ should be used in the fit to determine the coefficients a and 
b. For instance, the extracted values of the phase shift in the vicinity of the resonance 
{p/m^y ^ 0.8 neatly follow the straight line with the parameters a, b, which were determined 
from Eqs. ( I33l) . using the central values of m^, F^ in Eq. ( I32l) . 



Now, we are in a position to compare our results to those of Refs. J27|, |28|. The difference 



in the real part of the resonance pole position is small - both results agree with an accuracy 
of better than one per cent. The effect is larger in the imaginary part. However, one 
should keep in mind that the magnitude of the imaginary part is approximately 50 times 
smaller than the real part. As one concludes from Fig. |6l a relatively large effect in the 
imaginary part could be related, e.g., to the fact that the effective range plot js^not exactly 

28j and in 



linear. Therefore, it seems plausible that the systematic errors both in Refs. [2j 

the present paper are underestimated. We expect that the results should agree within the 

errors. 

Comparing our method to Liischer's approach, we further note that, once the plateau in 
L sets in, the energy and the width within our method can be extracted at a single value 
of L. In contrast to this, Liischer's approach implies the study of the volume dependence of 
the energy levels. This difference can be related to the fact that our method uses additional 
input information from MC simulations. In particular, apart from the energy spectrum, 
the two-point function Cii(t) contains the information about the transition matrix elements 
encoded in the constants zi, 1^1. 

Last but not least, we have also checked that our method works in the non-interacting 
case as well. Setting g = Q and adjusting k^, k^ in the Lagrangian to keep the masses of 



and 7] the same as in the interacting case (see Refs. 27|, |28|), we have done the calculation of 
the function Cii(t) anew. The fitted width turns out to be two orders of magnitude smaller 
as compared to the interacting case. This obviously corresponds to a stable particle. 
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(cf. with Fig. 2 of that paper). The errors in the data are purely statistical. In addition, we have 
indicated the dimension of the operator basis Oj used to extract the spectrum. 



VII. CONCLUSIONS 
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IV 



In the present paper we have proposed a novel method to extract the resonance pole 
position on the lattice. The method is based on the universal representation of the two- 
point function Eq. (I2U]) . which is valid given the sole assumption that an isolated low- 
lying resonance is present in the system. The energy and the width of this resonance 
are determined from the fit of Eq. fl20|) to the lattice MC data. It remains to be seen 
whether such a universal representation can be derived in more complicated cases (e.g., 
for the multi-channel scattering, see Ref. [38|) as well. 

The proposed method provides an alternative to Liischer's approach to the resonances. 
In the latter, the volume dependence of the spectrum on the moderately large lattices 
is studied. The spectrum consists of the scattering states only - the resonance has 
already decayed. In our approach the two-point function is studied at finite times t 
(when the resonance is still "alive") and for the asymptotically large values of L. Note 
that the actual calculations do not seem to require extraordinarily large volumes. For 
example, in the 1-1-1 dimensional Ising model L = 48 was already sufficient. 

The above difference entails an important advantage of the method described in this 
paper: whereas in Liischer's approach the MC simulations should be performed at 
least at several volumes in order to extract the resonance, the measurement at one, 
albeit sufficiently large, lattice volume suffices in our method. 

In certain cases, the numerical accuracy of the method can be improved considerably, 
if the multi-exponential representation of the two-point function is used in the fit 
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instead of calculating this function directly through MC simulations at all values of 
t. The coefficients of the multi-exponential representation are obtained by solving the 
generalized eigenvalue problem and, in particular, encode the decay matrix elements. 

v) Recently, the excited meson and baryon spectra have been determined by several lattice 



collaborations using the generalized eigenvalue equation (see, e.g. Refs. |19|-|21| for the 
latest work on the subject). These calculations closely resemble the calculations in 
the 1+1 dimensional toy model, which were presented in this paper. In our opinion, 
it would be very interesting to apply the proposed method to the data and if possible 
try to locate the resonance pole(s). This can be done at no additional cost, since the 
results of already existing MC simulations would be used. 
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Appendix A: Continuum limit in the matrix elements 



The regular summation theorem [29[, which is used in order to perform the continuum 
limit in the sums over the discrete momentum eigenvalues, implies that the integrand is a 
continuous function in the momenta. However, the finite- volume matrix elements (O|0(O)|/3), 
which enter Eq. ([8]), contain Lixscher's zeta- function and can become singular. Here, for one 
particular example, we shall demonstrate how these singularities are lifted. 

The averaged quantities, for which the validity of the regular summation theorem will be 
checked, are defined in Eq. f lTT]) . From now on, without loss of generality, we shall work in 
the center-of-mass frame k = 0. We wish to demonstrate that 

lim A(wo,A,0)=A°°(wo,A,0), wo>Wmin, A>0. (Al) 

L— >oo 

More precisely, the difference between the both sides of the above equation vanishes faster 
that any negative power of L, as L — )■ cxo. 

Let cuo be in the elastic scattering region. Since L is large, characteristic momenta are 
small and non-relativistic quantum mechanics provides an adequate description of a problem 
under consideration. Let us consider two massive (distinguishable) particles in the CM 
frame. The state vector corresponding to the eigenvalue Efj is given by 

L Yl M"^^ 1^' -^) ' 1^' -^) = «I(q)4(-q)|o) , (A2) 



ld/2 



q 



where aj, i = 1,2 denote the creation operators for the particles 1 and 2, respectively, and 
the wave function //^(q) is normalized, according to 

(/5i/3) = ^Ei//'(q)i' = i- (A3) 



q 
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For simplicity, let us further assume that the interaction between the particles is described 
by a separable potential l^(p, k) = gv{p)v{k), where the function v{p) corresponds to a 
smooth cutoff at large momenta. Note that in the following we will never need the explicit 
form of this function. The T-matrix is given by 

where // denotes the reduced mass of the system. 

In the limit E ^ Ep the T-matrix has a pole. In the vicinity of the pole, it behaves as 

Tfn a- E) = ^(P)^^^) = ^(P)^'(q) f A5l 

^P'^' ^ g~i_i(^Ep)-r{Ep){E-Ep) + ... -r[E„){E-Ep) + ...- ^ ^ 

From this expression, we may read off the wave function corresponding to the eigenvalue Ep 



The normalization of this wave function was chosen so that //^(p) obeys Eq. ( lA3p . 

Take the composite field 0(0) = 0i (0)02(0), where the 0i(O), i = 1,2 denote the elemen- 
tary particle fields. In momentum space, 

0.(O,x) = -^$^e^^''a.(k) (A7) 

k 

The matrix element that enters the spectral function is given by 

q 

where //3(r) denotes the Fourier-transform of /^(q). The averaged spectral function is writ- 
ten in the following form 

A{ujo, A, 0) = -1 5^ e{uo + A/2 - Ep)e{Ep - uo + A/2) \fp{0)\' . (A9) 

Hence, in order to verify the applicability of the regular summation theorem in this case, it 
suffices to show that //3(0) is a regular function of ko/s- As anticipated, this function contains 
Liischer's zeta-function which is singular at k^a = (27rn)^/L^. However, the factor I^Ep), 
which enters the normalization, contains the zeta-function as well. It is easy to check that the 
singular factors in the numerator and the denominator cancel, and the regular summation 
theorem holds. 



Appendix B: Emergence of the second Riemann sheet in the infinite volume limit 

Let us consider^ the function Fl of the complex variable z 



We are indebted to Jiirg Gasser who indicated this example to us. 
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Note that this function resembles the propagator of an unstable particle in the 1+1- 
dimensional effective field theory |5|. Further, it can be shown that 

sign(Im (A/icot(A/zL))) = — sign(Im (z)) . (B2) 

According to this condition, the denominator can not vanish outside the real axis. Thus, 
the only singularities of Fl{z) are simple poles on the positive real axis. 

If Im (z) 7^ 0, in the limit L — )• oo we have cot{y/zL) — )• — isign(Im (z)) and, therefore, 

Foo{z) = - ^\ ^ . (B3) 

1 — z -\- e^J—z 



The difference \Fi^{^z^ — Foo(2;)| vanishes exponentially with L, if Im (z) ^ 0. Note that, 
unlike Fl{z), which is a meromorphic function, Foo{z) is analytic in the complex plane cut 
along the positive real axis. This is what is meant when we speak of "the poles merging into 
the cut". 

Moreover, the function F^oiz) has a couple of complex-conjugated poles on the second 
Riemann sheet. These poles are solutions of the equation 1 — z + e^^—z that gives z± = 
1 =F ie^ + O(e^). If e^ is small, these poles come close to the physical scattering region and 
influence F^{z) on the physical sheet. Since away from the real axis the difference between 
Foo(-z) and -Pl(^) vanishes exponentially at a large L, the effect of the poles on the second 
Riemann sheet is felt in Fl{z) as well. 
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